On computing the coalescence time density in an isolation-with-migration model with few samples.

where f(t |Θ) is the probability density of the coalescent time. Assuming free recombination between loci, the full likelihood is a product of each locus likelihood. In order to determine the likelihood we must determine the density f(t |Θ) for coalescent of two samples in the IM model. Wang and Hey (2009) find the time to the MRCA by explicitly integrating over all possible sample paths in the system. The purpose of this letter is to demonstrate that the time to the MRCA is easy to compute from a matrix exponential. Furthermore, the matrix exponential framework has the advantage that it generalizes to more than two samples. Below we will describe the IM model, then briefly describe the solution to computing the coalescence time density from Wang and Hey (2009), and finally present an approach that computes the density through matrix exponentials. The description of the model is taken from Wang and Hey (2009) and formulates the IM model as a continuous time Markov chain. Before time T and for two samples the system is in one of the following five states: S11: Both genes are in population 1; S22: Both genes are in population 2; S12: One gene is in population 1 and the other is in population 2; S1: The genes have coalesced and the single gene is in population 1; S2: The genes have coalesced and the single gene is in population 2. The rates between the states can be described by the